
function [p q] = SI_Trap_Upwinding(dt)
global U  h  g  p  q  nt dx h1 x p_global tmax L steps itermax ;

dx1=1/dx;
dx2 = 1/(dx*dx);
dt2= 1/dt;


%---- Construct the spatial Discretization Matrix------%
A = zeros(x,x);
A(1,1) = U*dx1;
A(1,x) = -U*dx1;
for k=2:x
    A(k,k) = U*dx1;
    A(k,k-1) =  -U*dx1;
end

B=zeros(x,x);

B(1,x)=h*dx2;
B(1,2)=h*dx2;
B(1,1)=-2*h*dx2;
for k=2:x-1
    B(k,k)= -2*h*dx2; 
    B(k,k-1)=h*dx2;
    B(k,k+1)=h*dx2;
end
B(x,1)=h*dx2;
B(x,x-1)=h*dx2;
B(x,x)=-2*h*dx2;

C = g*eye(x);
I=eye(x);
beta1=0.5;
D = (dt2*I + beta1*A);
p_n=p;
q_n=q;

for n=2:nt+1; 
      
       for iter =1:itermax
           % advance q explicitly by Forward Euler method
           if (iter==1)               
               if (n==2)                   
                 qex(1:x,1) = -dt*(C*p_n(1:x,1)+A*q_n(1:x,1)) + q_n(1:x,1);
               else 
                  qex(1:x,1) = -2*dt*(C*p_n(1:x,1)+A*q_n(1:x,1)) + q_nm1(1:x,1);
               end
           end
           % solve p implicitly Trapezoidal
           
            rhs = dt2*p_n(1:x) - (1-beta1)*(A*p_n(1:x,1)+B*q_n(1:x,1)) - beta1*B*qex(1:x,1);
            p_np1 =D\rhs;
            
            % Correct q implicitly Trapezoidal            
            rhs = dt2*q_n(1:x) - (1-beta1)*(C*p_n(1:x,1)+A*q_n(1:x,1)) - beta1*C*p_np1(1:x,1);      
            qex  =  D\rhs;
       end
       q_nm1 = q;
       q=q_n;
       q_n   = qex;
       p_n   = p_np1;       
       p = p_n;
       
   
    if rem(n,100)==0

        refreshdata(h1,'caller') % Evaluate p in the function workspace
        drawnow
    end

end
display('Completed Successfully');